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ABSTRACT 

Submillimetre images of transition discs are expected to reflect the distribution of the 
optically thin dust. Former observation of three transition discs LkHa 330, SR 21N, and 
HD 1353444B at submillimetre wavelengths revealed images which cannot be modelled 
by a simple axisymmetric disc. We show that a large-scale ant icy clonic vortex that 
develops where the viscosity has a large gradient (e.g., at the edge of the disc dead 
zone), might be accountable for these large-scale asymmetries. We modelled the long- 
term evolution of vortices being triggered by the Rossby wave instability. We found 
that a horseshoe-shaped (azimuthal wavenumber m = 1) large-scale vortex forms by 
coalescing of smaller vortices within 5 x 10^ yr, and can survive on the disc life-time 
5 X 10^ yr), depending on the magnitude of global viscosity and the thickness of the 
viscosity gradient. The two-dimensional grid-based global disc simulations with local 
isothermal approximation and compressible-gas model have been done by the GPU 
version of hydrodynamic code FARGO (GFARGO). To calculate the dust continuum 
image at submillimetre wavelengths, we combined our hydrodynamical results with a 
3D radiative transfer code. By the striking similarities of the calculated and observed 
submillimetre images, we suggest that the three transition discs can be modelled by 
a disc possessing a large-scale vortex formed near the disc dead zone edge. Since the 
larger dust grains (larger than mm in size) are collected in these vortices, the non- 
axisymmetric submillimetre images of the above transition discs might be interpreted 
as active planet and planetesimal forming regions situated far (> 50 AU) from the 
central stars. 

Key words: accretion, accretion discs — protoplanetary discs — hydrodynamics — 
instabilities — methods: numerical) 



1 INTRODUCTION 



Recently Brown et al. ( 2009 ) presented submillimetre images 
of three transition discs LkHa330, SR21N, HD 135344B. 
Using a 2D radiative transfer model (|Dullemond & Dominik] 
2004 ) , they showed that the observed inner discs are cleared 
from dust as the optical through millimetre-wave spectral 
energy distribution (SED) can be fitted with a disc having an 
inner cavity. The non-axisymmetric intensity distributions 
(see fig. 1 of Brown et al. 2009), however, require further 



assumptions such as the presence of an embedded planet at 
large distances from the central star (> 50 AU) that suffi- 
ciently perturbs the disc dust distribution (see, e.g., Wyatt 



et al. (1999)). A massive enough planet opens a gap, and 



the disc's material will be piled up just outside of the gi- 
ant planet's orbit ( de Val-Borro et aL]|2007 ). Although the 
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giant planet induced non-axisymmetric density distribution 
hypothesis is challenging, we propose in this Paper that the 
observed asymmetries are due to large observable vortices 
formed in these discs during their early evolutionary stages. 

However, the current most accepted scenarios for giant 
planet formation, the core-accretion (Bo denheimer Pol-] 
lack||1986 Pollack et al.||1996 ) and the gravitational insta- 
bility ( Boss||20^ 2006) face problems forming and keeping 
giant planets at such a large distances. Within the core- 
accretion scenario, the optimal places for planet formation 
is the vicinity of the snow line where the condensation of 
volatiles into ices increases the surface density of solids by 
a large factor, increasing the isolation mass above the criti- 
cal core mass. The location of the snow line determined by 
the opacity of the dust grains and the mass accretion rate 
is recently found to be < 16 AU ( |Min et al.||2011t . Eor sig- 
nificantly larger distances, the growth rate of the planetary 
core is strongly limited. On the other hand, forming giant 



2 Zs. Regdly, A. Juhdsz, Zs. Sdndor and C. P. Dullemond 



planets on wide orbits via gravitational instability does not 
face many problems (assuming rapid cooling of the gas), 
except that it is not yet clear how to keep them there. 
[Dodson-Robinson et al. (2009) investigated three possible 
mechanisms - core-accretion (with or without migration), 
scattering from the inner disc, or gravitational instability - 
which might result in massive planets at wide (> 35 AU) or- 
bits discovered by high-contrast direct imaging observations 
of Fo malhaut (K alas k et al.,|2008J, HR8799 (Marois et al.^ 



Kretke & Lin ( 2007 ) have already shown that pressure- 



2008) or CD-352722 (Wahhaj & etal., 2011) for instance. 



Dodson-Robinson et al. found that the only possible sce- 
nario to form stable wide-orbit planets is the gravitational 
instability, but planets formed this way were recently shown 
(by Baruteau et al.|20TT ) to migrate quite rapidly inwards - 
puts in question its feasibility. Note that ( Crida et al.|2009 ) 
have shown that runaway outward migration by planets in 
mean motion resonance carving a common gap in the proto- 
planetary disc is an other possible mechanism to place giant 
planets on wide orbits, if the inner one is significantly more 
massive than the outer one. 

In recent years the formation of vortices in protoplane- 
tary discs have been investigated thoroughly, since they may 
play an important role in the formation of planets in the 
core-accretion scenario. Observation of such vortices would 
have a great importance providing further informations and 



support for planetary formation theories (see, e.g., Klahr & 



Henning| ( p97J ; Barge Sommeria (1995) ). 



The classical approach of accretion discs (Shakura & 
Sunyaev]|1973 ) postulates a simple scaling relation between 



the gas viscosity and the level of angular momentum trans- 
port parametrized by a. An accretion disc coupled dynam- 
ically to a weak magnetic field is subject to the so-called 
magnetorotational instability (MRI) resulting in outward 
angular momentum transport, which can be described by 



an effective a-viscosity (Balbus & Hawley 1991|. The disc 



becomes MRI-active, if the ionisation state of the gas is suf- 
ficient. Investigating the ionisation processes in MRI-active 
discs, Gammie ( 1996 ) proposed the existence of a radially 



extended region in the disc midplane, where the accretion 
is strongly reduced. Near the disc inner edge, the midplane 
gas is collisionally ionised, owing to the high temperatures 
( Pneuman Mitchell|1965l|Umebayashi|1983| ). Further out, 
the temperature is not high enough to collisionally ionise the 
gas, while the gas is yet dense enough to shield the interstel- 
lar cosmic-rays or diluted stellar X-rays. The magnetic field 
of the ionised disc corona might even deflect the cosmic rays. 
As a result, only a thin surface layer of the disc is ionised 
at this region, while at even larger distances, the cosmic 
rays can penetrate the tenuous disc which might ionise the 
disc midplane again ( Glassgold et al.|1997 Igea & Glassgold 



1999 ). As a consequence, the turbulent viscosity is low in the 



poorly ionised region of the disc midplane, resulting in the 
formation of the so called dead zone. Within this region the 
viscous accretion is taking place mostly in a narrow upper 
layer of the disc, while it is stalled in the disc's midplane. 
Varniere k Tagger] ( |2006| ) and |Terquem| ( |2008| ) demon- 



strated that density and pressure enhancements ('bumps'^ 
appear at the boundaries of a dead zone because of the sharp 
jump in the viscosity. These pressure bumps are found to 



be unstable to the Rossby wave instability (RWI) ( Lovelace 



et al. 1999 Li et al. 2000), leading to the formation of a 



gradient inversion formed at the inner edge of the dead zone 
or at the snow line causes the accumulation of dust. The 
existence of large-scale vortices can also accelerate the for- 
mation of planetesimals and planetary embryos in the core- 



accretion paradigm. Barge & Sommeria ( 1995 ) and Klahr & 



Henning (1997) have already shown that particles tend to 



get trapped in anticyclonic vortices. At the radial location 
of the vortex the type-I migration is halted due to radial 
pressure-gradient inversion, thus not only dust is retained 
here but larger bodies too. Particles of about centimetre 
to metre in size subsequently drift into the pressure max- 
ima, i.e., to the vor tices |Lyra et al.|2 009; Kato et al.|2010l 
Dzyurkevich et al. 2010| and form gravitationally bound 



clumps of solids, which can coalesce forming embryos be- 
tween the masses of the Moon and Mars. The swarm of these 
embryos evolves further by mutual N-body collisions form- 
ing massive cores (~ IOM0) of giant planets in ^ 5 x 10^ yr 
(Sand or et al.||2011| . Thus, large-scale anticyclonic vortices 
formed in discs are extremely important for the planetary 
formation helping to overcome the meter-size barrier prob- 
lem (Blum & Wurm 2008), and the time-scale problem of 



oligarchic growth. From this perspective these vortices can 
be regarded as 'planetary factories'. 



The micron sized dust is dynamically coupled to the 
gas as it is indirectly supported against the settling by gas 
pressure via the gas drag, thus their spatial distribution is 
expected to be the same. Assuming that the three transi- 



tion discs presented in Brown et al. (2009) contain signif- 
icant amount of gas - supported by Ha (Fernandez et al. 



1995) , [01] (van der Plas et al.|2008t and fundamental ba nd 
CO emissionT p^ontoppidan et al.|2008||Salyk et al.|2009| ) -, 

these vortices can be imaged in submillimetre wavelengths 
as overdense regions in the dust. The inner holes in the sub- 
millimetre images of the three discs may indicate that dust 
is either consumed here by a planet forming process or by 
coagulation to larger sized grains, which are drifted rapidly 
into the star. Note that coagulation without fragmentation 
would make the vortex also invisible as the dust opacity is 
decreasing with grain size. If the dust trapped inside giant 



large-scale anticyclonic vortices. 



vortices (Inaba Barge 2006) is accountable for the ob- 
served features, fragmentation should be happening inside 
them, in order to replenish the population of small dust 
grains. 



Inspired by these ideas, we explore, whether a large- 
scale, long-lived (on 10^ yr time-scale) non-axisymmetric 
gas-dust perturbation formed at the outer edge of the disc 
dead zone could explain the non-axisymmetric submillime- 
tre images. To calculate the density perturbation, we run 
long-term numerical hydrodynamical simulations on a GPU 
aided computer. In § 2, we present our hydrodynamical disc 
model, and a study on the robustness of the formation of 
large-scale vortices. Submillimetre dust continuum images 
calculated by a 3D radiative transfer using our hydrody- 
namical results are presented in §3. The Paper closes (§4) 
with the discussion and concluding remarks, along with a 
short discussion on the reasonability of our findings. 
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Figure 1. Example of smooth viscosity reduction according to 
Eq. ([T]) used in our models. Sa{R) for models a = 0.01, amod = 
0.01, Rdze = 50 AU, and Ai^^ze = 1? 3, 5 AU are shown with black 
solid, blue dashed, and red dotted curves, respectively. 



2 VORTEX FORMATION 



We revisit the development of RWI, ( Lovelace et al.|p"999 ), 
by 2D, grid-based, global hydro dynamical disc simulations, 



assuming the a-type viscosity prescription ( Shakura & Sun- 
To model the dead zone 



yaev 



1973) 



smoothly reduced such that a{R) 



the value of a, is 
aSa{R) and 



- 0.5 (1 - ttmod) 



- tanh 



( R — Rdze \ 
V Ai?dze J 



+ 1 



(1) 



where ttmod is the depth of turbulent viscosity reduction, 
Rdze and ARdze are the radial distance and the width of 
the transition region of turbulent viscosity at the dead zone 
outer edge, respectively. Regarding the present study the 
effect of the inner edge of the dead zone is not relevant, 
therefore it is not investigated here. Examples of a{R) used 
in our models are shown in Fig.^ 

To quantify i?dze, we adopt the results of |Matsumura| 
& Pudritz (2005). According to their calculations, the outer 



edge of the dead zone lies between 12AU and 36 AU, de- 
pending on their various model parameters. To induce the 
vortex formation at distances suggested by submillimetre 
images of ( [Brown et aLl[2009t , we set the outer edge of the 
dead zone to slightly larger distance (i^dze = 50 AU). Note 
that the exact location of the outer edge of the dead zone 
is very uncertain, because the physics of the ionisation (de- 
pending on the ionisation rate of the stellar X-rays, interstel- 
lar cosmic-rays, proton radiation from the high temperature 
disc corona, decay of radionuclides such as ^"^K or ^^Al) and 
the dependence of the magnetorotational instability on this 



ionisation degree is not yet well understood (see, e.g.. Turner 
k Drake||2Q09t . 



2.1 Numerical simulations 



To follow the evolution of the RWI triggered vortices on 
several 10^ yr time-scale, we used GFA RGCQ the GPlQver- 
sion of the FARGO code ( |Masset|20QQt . We slightly modified 
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^ Graphical Processor Unit constructed with several hundred 
procesors (240 for GTX-280, which we used for the simulations). 
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Figure 2. The distribution of Qh (where Q is the Toomre pa- 
rameter and h is the disc aspect ratio) calculated for the base 
model, when a large-scale vortex is fully-developed. Since Qh is 
well above unity, even inside the vortex eye, seems reasonable to 
neglect the self-gravity. 



the CUDAE] kernel of GFARGO code in order to handle the 
viscosity reduction according to Eq. ([T]) . 

The dimensionless continuity and Navier-Stokes equa- 
tions , describing the surface mass density evolution in an 
accretion disc, are numerically solved in two-dimensional 
cylindrical coordinate system. We use the locally isother- 
mal approximation, in which the vertically integrated gas 
pressure is p{R) = Cs(i?)^S(i?), where Cs{R) is the local 
sound speed, and S(i?) is the vertically integrated surface 
mass density of gas. The initial surface mass density pro- 
file follows the power- law distribution, S(i?) - R'^-^, We 
use radially constant disc aspect ratio h = H{R)/R, where 
H(R) is the pressure scale height at the radial distance i?, 
for which case Cs{R) = QH = y^GM./R^{hR) oc R'^^^. 

The disc self-gravity can be neglected as long as 
the Toomre parameter ( Toomre] [l964) is high, Q = 
/iM*/(7ri?^S) > 1, which holds in our unperturbed discs. 
Note, however, that Q is relevant only for small scale, i.e. 
high m modes, for global modes (e.g., m=l) the parameter 
that might really measure the relevance of self-gravity is Qh 
( fPapaloizou 2002). Since Qh might be smaller than unity. 



particularly inside the vortex owing to the density enhance- 
ment, we calculate Qh for the base model {h = 0.05 and 
Mdisc = 0.01 M0, and neglecting self-gravity) by the time 
when the large-scale vortex have fully developed, see Fig. [2] 
As one can see, Qh never reaches unity (Q/i|min = 3.68 inside 
the vortex eye), therefore we neglected the disc self-gravity 
in our models. Note that the disc mass Mdisc = 0.01 Mq that 
we used in our simulations is in the range (0.006 ^ Mdisc ^ 
0.026) as has been derived by SED modelling for the three 
transition discs by [Andrews et al. (2011). Nevertheless, Qh 
might be smaller than unity for a disc which is three times 
more massive than the one we used, thus for that case one 
should take into account the self-gravity. 

The disc inner and outer boundaries are fixed at 4AU 
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and 100 AU, respectively. The computational domain con- 
tains 256 radial logarithmically and 512 azimuthal uniformly 
distributed grid cells. The numerical convergence with the 
above resolution is found to be satisfactory, see Sect. |2.3| 
At the inner and outer boundaries the so called 'damp- 



100 



ing boundary condition' are applied (de Val-Borro & etal., 



2006 ) , where a solid wall is assumed with wave killing zones 



next to the boundaries. 



2.2 Development of a large-scale vortex 

In what follows we give a short description of the formation 
and evolution of a large-scale vortex in our base disc model 
(Mdisc = 0.01 Mo, h = 0.05, a = 0.01, a^od = 0.001, i?dze = 
50 AU, and ARdze = lAU). The major evolutionary stages 
of the surface mass density distribution is shown in Fig|3] 
while the azimuthally averaged surface mass density profiles 
are shown in Fig.|4] 

To excite RWI, a small mass orbiting body (negligible 
to the local disc mass) was put in the disc at 50 AU. The or- 
biting perturber excites low-m mode density waves emerging 
from its Lindblad resonances, and these in turn trigger the 
RWI. Triggering RWI by this method is appropriate, because 
RWI is a low-m mode instability. We observed, though, that 
for models where the RWI is strong - e.g., for sharp dead 
zone edge models (Ai?dze = lAU), where the radial gradi- 
ents of pressure and density are steep - the numerical noise 
in the radial velocity component can be large enough to 
trigger the vortex formation. This can be explained by the 
fact that the magnitude of radial velocity component is very 
small in Keplerian discs, thus a small numerical noise (which 
is in the order of 10~^ in the radial velocity component for 
our simulations done by GFARGO) is enough to excite RWI. 
For this reason, we investigated whether the mode of trig- 
gering RWI (by perturbing the radial velocity component or 
the surface mass density with superposition of sines or white 
noise; igniting density waves by a small mass planet orbiting 
at 50 AU) has any effect on the vortex evolution. Studying 
carefully the results of our simulations, we have not found 
any significant differences, see details in Appendix B. 

According to our simulation, the gas accumulates soon 
(^ 10^ yr) at slightly smaller distance than the outer edge 
of the dead zone (see, e.g.. Fig. [4] solid line curve). The RWI 
is excited at 4 x 10^ yr (12 orbits at 50 AU assuming a 
Solar mass central star, i.e. Pso,! = 12) forming several 
anticyclonic vortices with azimuthal wavenumber m = 6 
(Fig.js] Panel A). These vortices coalesce, while the density 
bump continuously strengthens (see, e.g.. Fig. [4] dashed line 
curves). As a result, only one large-scale vortex (m = 1 mode 
vortex) remains by ^ 1.5 x lO'^ yr (Pso,! = 40) (Fig.|3] Panel 
B). Both the azimuthal and radial extensions of vortex in- 
crease, producing a horseshoe-shaped density perturbation 
by - 10^ yr (Pso,! = 300) (Fig|3]panel C). The vortex wings 
finally locks by - 3 x 10^ yr (Pso,! = 900) (Fig.^ Panel 

D) . Investigating the following evolution of large-scale vor- 
tex in the base model, we found that the horseshoe-shaped 
density perturbation disappears, and a ring like density per- 
turbation develops at 5 x 10^ yr (Pso,! = 1400) (Fig.|3] Panel 

E) . On the top of the ring- like density perturbation, small- 
scale vortices form, which are also subject to the vortex co- 
alescing process (see in the next paragraph). After ^ 10^ yr 
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Figure 3. Major evolutionary stages of the surface mass density 
evolution in our base model {MfUsc = 0.01 Mq, H/R = 0.05, a = 
0.01, Ai^dze = 1 AU, amod = 0.001, ). Panel A: excitation of the 
RWI at ~ 4 X lO^yr (Pso,! = 12). Panel B: only one large-scale 
vortex remains after ~ 1.5 x 10^ yr (Pso,! = 40). Panel C: vortex 
evolution results in the formation of horseshoe-shaped (azimuthal 
wave number m = 1 mode) surface mass density distribution at 
~ 10^ yr (Pso,! = 300). Panel D: vortex wings are interlocked 
at ~ 3 X 10^ yr (Pso,! = 900). Panel E : ring-like structure forms 
with small vortices on the top of that at ~ 5 x 10^ yr (Pso,! = 
1400). Panel F: reappearance of the horseshoe-shaped vortex at 
-6x lO^yr (Pso,! = 1700). 



(Pso,! = 300) of evolution, the m — 1 mode (horseshoe- 
shaped) vortex reappears (Fig.[3] Panel F). 

To illustrate this vortex mode oscillation we calculated 
the azimuthal spectral power of the surface mass density 
distribution against time measured in the orbital period 
Pso,! shown in Fig[5] First, a periodic Fourier transform of 
the density in the ring (period of 27r along ^-direction at 
the radius where the vortices are) was calculated. Then the 
Fourier amplitudes of the m = — 6 modes are calculated 
normalised by the total power in these seven modes. It is 
clearly seen that the spectral power is mainly in the m — 1 
mode (black curve) throughout the whole simulation. How- 
ever, by 5 X 10^ yr (Pso,! = 1400) of disc evolution, the 
spectral power in m = 1 mode temporarily drops. When the 
m = 1 mode vortex decays, the spectral power of m = 4, 5 
modes abruptly rises, i.e. RWI is re-excited. These newly 
formed vortices are also subject to coalescing and form again 
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Figure 4. Evolution of the azimuthally averaged surface mass 
density profile in the base model. The solid black curve represents 
the initial surface mass density profile. The density profiles shown 
with coloured dashed curves are calculated from the surface mass 
density distribution presented in Fig[3] 
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Figure 5. Azimuthal spectral power of the surface mass density 
distribution against time measured in Pso,! for model presented in 
Figjs] The Fourier amplitudes are normalised by the total power 
measured in the m = — 6 modes. 



an m = 1 mode vortex by - 10^ yr (Pso,! = 300). Thus, 
the typical vortex mode is m = 1 through 1.8 x 10^ yr 
(^50,1 = 5200) of vortex evolution, which is interrupted by 
short-lived higher mode (1 < m < 5) vortex configurations. 
Note that the spectral power in m = characterising the 
average surface mass density drops off immediately as RWI 
is excited. 

As shown in Fig.^ the density maximum is pushed in- 
ward from its original location, i.e., the m — 1 vortex is 
slightly shrinked due to the extra pressure caused by the 
density build-up. As a result, the fully-developed m — 1 
mode vortex forms at slightly smaller distance {R 40 AU) 
to the central star than the distance where the viscosity re- 
duction (Pdze = 50 AU) begins. It is also appreciable that 
the density enhancement inside the vortex eye increases 
with time. Finally, emphasise that the surface mass den- 
sity of the disc is significantly lowered next to the dead 
zone edge. This density rarefaction continuously strengthens 
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Figure 6. Fully-developed large-scale vortices at 3.5 x 10^ yr 
(Pso,! = 1000) in models using Ai^^ze = 1? 2, and 3 AU dead zone 
edge width shown in panels A, B, and C, respectively. The ax- 
isymmetric m = configuration for model Ai^^ze = 5 AU shown 
in panel D is just the result of the viscosity reduction at the dead 
zone edge stalling the mass infiow, with subsequent pileup of gas. 



with time, forming a deep hole in the disc by ^ 1.5 x 10^ yr 
(Pso,! = 4250). 



2.3 Robustness of large-scale vortex development 

To examine the robustness (sensitivity to disc properties) 
of the development of the large-scale vortex, we performed 
several hydrodynamic simulations. By harnessing the power 
of the GPU, each simulation required only a day of computer 
time. For this parameter study the boundary conditions, the 
disc mass, and the disc aspect ratio were the same as in the 
previously shown base model. In what follows we summarise 
briefly the results of these simulations. 



2.3.1 Numerical resolution 

To check whether our numerical solution is satisfactory with 
the original grid {Nr — 256, Ncp — 512), we recalculated the 
base model with higher resolution, namely Nr = 512, A/'^ = 
1024. We found that the results are practically identical, 
even by frame-to-frame comparison. Hence, the spatial res- 
olution of the base model (256x512) is sufficient for the 
following simulations. 



2.3.2 Numerical precision 

Scientific calculations with conventional GPU cards (e.g., 
NVIDIA GTX-280 that we used) lacks the usual double pre- 
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Figure 7. Azimuthally averaged surface mass density profiles at 
3.5 X 10^ yr (Pso,! = 1000) in models using Ai? = 1, 2, 3, and 5 AU 
dead zone edge width and OLmod = 0.001 {upper panel). Models 
with lower viscosity reduction amod = 0.01 (lower panel) show 
no significant departures to previous simulations. 



cision at the moment |^ To check the validity of our results 
provided by the single precision GPU we compare them to 
those of given by the CPU based FARGO code in the fastest 
evolving model (a = 0.1, ARdze = 1 AU, ttmod = 0.001). 
The two runs give qualitatively the same result: the evolu- 
tion of the large-scale vortex (growth rate, vortex coalescing, 
wing interlock, and oscillation) is found to be practically the 
same (see the details in Appendix A). 



^ 10"^ 



< m = 

m = l 

m = 2 - - 

\\ m = 3 -■- 

\ \ , m = 4 

\'\ \^ m = 5 

■1 ,m = 6 — 



, \ / \ h\ ;/^ i; i ;\ 



,'1 If i 



1 



2.0x10^ 4.0x10^ 6.0x10^ 8.0x10^ I.OxlO'* 1.2x10^ 1.4x10^ 
Time (P50,,) 

Figure 8. Azimuthal spectral power of the surface mass density 
distribution against time measured in Pso,i for model a = 0.01, 
Ai^dze = 3AU, 6a = 0.01. The Fourirer amplitudes are nor- 
malised by the total power measured in the m = — 6 modes. 



stalling the mass inflow, with subsequent pileup of gas. Thus, 
the RWI excitation occurs only in ARdze ^ 2i7 models in a 
good agreement with [Lyra et al. (2009). 

Figures|6] and [7] show the surface mass density contours 
and the azimuthally averaged surface mass density profiles 
developed at 3.5x10^ yr (Pgo,! = 1000) for various dead zone 
edge width models. Although the vortex structures are dif- 
ferent (the larger the dead zone edge width, the smaller the 
azimuthal extension of the m — 1 mode vortex, see Fig. [6] 
Panel A, B, C), the azimuthally averaged surface mass den- 
sity profiles are very similar (Fig.[7| dashed line curves). 

The vortex mode oscillation behaviour is completely ab- 
sent in models with ARdze ^ 2 AU during the whole simu- 
lation covering 5 x 10^ yr (Pso,! = 14000). Figurejs] show 
the azimuthal spectral power of the surface mass density 
distribution against time calculated as described before. For 
this case, the largest spectral power is in m = 1 mode and 
remains nearly constant throughout the whole simulation 
after the initial vortex coalescing process. 



2.3.3 Dead zone 



width 



Li et al. (2000) have found that the excitation of the RWI 



occurs later in time if the transition region of the viscosity 
reduction (i.e., the width of the dead zone edge) is wider. We 
confirmed this by setting the dead zone edge width ARdze 
between 1 AU — 5 AU for models using a = 0.01 and ttmod = 
0.01. The excitation of the RWI occurs about 4, 8.5, and 
20 X 10^ yr (Pso,! = 12, 24 and 57) for models with 1, 2, and 
3 AU dead zone edge width, respectively. However, in model 
with Ai^dze = 5AU, the RWI is failed to be excited and 
only axisymmetric (m = 0) density perturbation develops 
as a result of the viscosity reduction at the dead zone edge 



^ Note that there exist scientific grade graphics cards in the 
NVIDIA Tesla series which are capable of double precision calcu- 
lations. 



2.3.4 Depth of viscosity reduction 

We investigated the effect of the magnitude of amod being in 
the range of 0.001 — 0.01 for models with a = 0.01. We found 
that the time required to ignite the RWI, and to form the 
fully-developed large-scale vortex are independent of amod- 
Moreover, the magnitude of density contrast developed in 
the large-scale vortex is also independent of amod- As a 
consequence, the azimuthally averaged surface mass density 
profiles (Fig.|7| calculated at - 3.5 x 10^ yr (Pso,! = 1000) in 
weak viscosity reduction (ttmod = 0.01, dashed line curves)) 
models are very similar to that of in models with stronger 
viscosity reduction (ttmod = 0.001, solid line curves)). The 
life-time of the m — 1 vortex, i.e., the vortex mode oscil- 
lation behaviour, is also found to be independent of amod^ 
for sharp dead zone edge models the vortex mode oscillation 
do occur, while for thick dead zone edge models the m — 1 
mode vortex is long-lived. 
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2.3.5 Magnitude of global viscosity 

We found that the vortex evolution is sensitive to the magni- 
tude of the global viscosity a being in the range of 0.001—0.1. 
The time required to excite RWI and to form the fully- 
developed large-scale m = 1 mode vortex decreases with 
increasing a. Both the radial extension of the large-scale 
vortex and the density contrast in the vortex increase with 
a. These can be interpreted by the accelerated vortex evo- 
lution in models with high global viscosity. 

For high viscosity models [a — 0.1), the inner disc be- 
comes ever denser during the later (t ^ 5 x 10^ yr) epochs, 
and the density contrast drops between the large-scale vor- 
tex and the inner disc. In consequence, the m — 1 mode 
vortex survives only up to ^ 10^ yr in these models. For low 
viscosity discs (a = 0.001), the vortex evolution slows down 
considerably, thus the excitation of RWI occurs later and 
the life-time of m = 1 mode vortex can reach ^ 5 x 10^ yr. 

We found that for low viscosity discs the vortex mode 
oscillation is completely missing, even for sharp dead zone 
edge (Ai?dze = 1 AU) models. Contrary to this, for high vis- 
cosity discs, the vortex mode oscillation does occur even for 
thick dead zones ( Ai?dze < 5 AU) models. Emphasise that 
the m — 1 mode (horseshoe-shaped) vortex might survive 
up to the life-time of the disc for a = 0.001, independently 
of the dead zone edge width being smaller than the twice of 
the local disc height (Ai?dze ^ 5AU). 



3 SUBMILLIMETRE IMAGES 

In order to show the resemblance of the observations of 
Brown et al. (|2009) and signatures of large-scale Rossby- 
vortices in protoplanetary discs, we calculated submillimetre 
images using the 3D radiative transfer code RADMC-3E[^ 
First, the temperature structure of the disc was determined 
in a Monte Carlo simulation then the images at 340 GHz 
(880 /xm ) were calculated via raytracing. Submillimetre Ar- 
ray (SMA) observations were simulated in the Common As- 
tronomy Software Applications package (CASA)[^ 

Snapshots of surface mass density distributions shown 
in Fig. |9] from the hydrodynamic simulations were converted 
into volume density assuming a vertically Gaussian density 
distribution with the same constant aspect ratio {h — 0.05) 
as used in the hydrodynamic simulations. The geometrical 
dimensions of the sources as well as their distances were 
adapted from Andrews et al.| (|2011| . The hydrodynamical 
frames had to be rotated and spatially scaled to approxi- 
mately match the disc shapes and sizes on the measured 
submillimetre images. The used rotation angles were 330°, 
320°, and 170° clockwise, while the linear scaling factors 
were 2, 1.2, and 1.7 for LkHa330, SR21N, and HD 135344B, 
respectively. Since there are evidences for a fairly large gap 
in all discs ( Brown et al.|2009" ) , we reduced the dust density 
within a certain radius according to Andrews et al. (2011). 



Inside the dead zone, dust grains may grow and drift rapidly 
inwards leaving a reduced density region, i.e. a gap, behind. 
In ant icy clonic vortices, however, dust grains get trapped 



^ http:/ /www. it a. uni-heidelberg.de/~dullemond/software/radmc- 
3d 

^ http:/ /casa.nrao.edu/index.shtml 



and cannot drift inwards, which would qualitatively explain 
why the outer radius of the gap is close to the location of 
the RWI. Note that there is also evidence for the presence 
of an inner disc in th ese sources (Fernande^ et al.|1995| [van 



der Plas et"ar][2008l [Pontoppidan et al.||2008| |Salyk et aT 



2009| ). Since these inner discs have an outer radius probably 



smaller than the inner radius of our hydrodynamic simula- 
tions (4 AU) we neglected them. 

In the radiative transfer simulations we used the hy- 
drodynamic grid in the radial and azimuthal directions 
while N6» = 100 grid points were used in the direction. 
The polar grid was a uniformly distributed, but we split 
the — 7r/2 ^ ^ ^ +7r/2 domain into two regions in order 
to resolve the disc vertically, while keeping the simulation 
computationally feasible. We placed 80 grid points in the 
1^1 ^ 0.35 domain and 20 grid points were placed in the 
1^1 ^ 0.35 region. In the thermal Monte Carlo simulations 
2 X 10^ photons were used to determine the temperature 
structure. 

The dust in our model consisted of "astronomical" sili- 
cate grains (Weinga rtner k, Draine|200T ), whose absorption 
and scattering cross sections were calculated using Mie the- 
ory from the optical constants. D ust grains had an MRN size 
distribution (Mathis et al. 1977), n(a) oc a~^"^, sampled with 



10 logarithmic spaced bins with minimum and maximum 
grain size of 0.1 /im and 1000 /im, respectively. Dust par- 
ticles with different sizes have different temperatures, thus 
they are not thermally coupled to each other. The dust-to- 
gas ratio was assumed to be 0.01 and constant within the 
disc outwards of the outer radius of the gap. Inside the gap 
the dust-to-gas ratio should be lower than 0.01 if the gap is 
formed due to the growth and inward drift of dust grains. 
Finally, the simulation of the observation (Fourier transfor- 
mation, u-v sampling, cleaning) was done with CASA. The 
resulting images are shown in Fig.[lO] We also calculated 
SED from our synthetic disc model for all three sources and 



compared them to the observations (see Fig. 11 ) 



Our synthetic model images (Fig. 10) resemble the most 



important features of the submillimetre continuum map of 
the three transition discs presented in Brown et al. (2009). 



The disc emission is clearly non-axisymmetric, horseshoe- 
shaped, and it shows single arc- like maximum at position 
angle around 225° for LkHa 330, 120° for HD 135344B, while 
two maxima are visible at 120° and 270° position angles for 



SR21N. The SEDs (Fig.[TT]) calculated from our disc mod- 
els match the observations shortwards of about 1 /xm and 
longwards of about 10 /xm. Between 1 and 10 /xm, however, 
our models under predict the observed fluxes, which is a 
natural consequence of neglecting the inner disc emission. 
We concluded that our disc models with large-scale Rossby 
vortices fit the observations of the mid-infrared to submil- 
limetre emitting regions of all three discs. 

Given the fact that our model is vertically optically thin 
at submillimetre wavelengths the observed asymmetry in the 
continuum emission can arise from temperature, density per- 
turbations or the combination of both. Since our disc is op- 
tically thin at 880 /xm the increase in the column density 
will increase the observed flux. The higher density regions, 
however, can be shielded from the heating radiation (higher 
optical depth) and therefore the temperature can decrease. 
The temperature perturbations in our model are less than 
a factor of two, which means the same factor for the flux 
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Figure 9. Hydrodynamic frames that were used to calculate the 
synthetic submillimetre images presented in Fig.^O] For mod- 
elling LkHo; 330 and SR21 frame shown in upper panels while 
for HD 135344B frame shown in lower panel were used. The 
frames were calculated in two different models: parameters were 
i^dze = 50 AU, Ai^dze = 1 AU, a = 0.01, a^^od = 0.001 for upper 
panel, and i?dze = 50 AU, Ai^dze = 2AU, o; = 0.1, a^od = 0.01) 
for lower panel. To match the observed submillimetre images of 
[Brown et al. (2009) frames had to rotate by 330°, 320°, and 
170° clockwise and spatially scale by 2, 1.2 and 1.7 for LkHa 330, 
SR21N, and HD 135344B, respectively. 



perturbation caused by temperature effects, since we are in 
the Ray leigh- Jeans regime. The density perturbations on the 
other hand are much larger in amplitude (factor of 6-8, see, 
e.g., the red dashed curve in Fig.[3]) and are therefore most 
likely responsible for the observed asymmetries in the im- 
age. A detailed study of these effects in non-axisymmetric 
discs will be the topic of a forthcoming paper. 



4 DISCUSSION & CONCLUSION 

In this Paper we propose an interesting physical phe- 
nomenon that can be responsible for the large non- 
axisymmetric density perturbations revealed by submil- 
limetre images of three transition discs of [Brown et al.| 
(2009). Based on the results of our long-term, two- 
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Figure 10. Submilhmetre images of LkHa 330 (top), SR21N 
(middle) and HD 135344B (bottom) calculated by the 3D radia- 
tive transfer code, from the frames of the hydrodynamical calcu- 
lations presented in Fig.|9] The synthesised beam is shown in the 
lower left corners in blue. The position angle of the discs, i.e. the 
vortex eye position was rotated corresponding to the observation 
of [Brown et al.| ( |2009| >. 
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Figure 11. The observed (black dots) and the theoretical SEDs calculated in the ray-tracing models of LkHaSSO (left), SR21N (middle), 
and HD 135344B (right). The failure to fit the NIR fiux is due to the fact that, for computational feasibility reasons, the inner disc emission 
is neglected. 



dimensional, grid-based, global hydrodynamic simulations 
(using compressible-gas model, neglecting the thermody- 
namical effects, and the disc self-gravity), we suggest that 
the observed asymmetries are caused by a large-scale vortex 
developed where the turbulent viscosity has a steep gradi- 
ent, e.g. at the outer edge of the disc dead zone. Vortices are 
triggered by the Rossby wave instability (RWI), which later 
are subject to a coalescing process and form a large-scale 
anticyclonic vortex. Note that [Wolf Klahr (2002) have 



previously shown that such large-scale anticyclonic vortices 
formed in the inner disc (< 10 AU) can causes detectable 
intensity asymmetry in millimetre wavelengths (345 GHz 
and 900 GHz) which are proposed to be observed by the 
fully operable Atacama Large Millimeter/submillimeter Ar- 
ray (ALMA) with the longest available baseline. 

After the initial excitation of RWI, a large-scale anticy- 
clonic vortex develops by the coalescing of small-scale vor- 
tices at ^ 50 AU, near the outer edge of disc dead zone, 
(Fig.|3]and Fig.|4|, which grows both in size and density con- 
trast. The time required to form a m = 1 mode large-scale 
vortex is found to be ^ 1.5 - 4 x 10^ yr (Pso,! = 40 - 115), 
depending on the dead zone edge width ARdze (the sharper 
the dead zone edge, the faster the vortex evolution) and the 
magnitude of global viscosity (the larger the viscosity, the 
faster the vortex evolution). The m = 1 mode vortex tem- 
porarily disappears and reforms through RWI re-excitation 
and subsequent vortex coalescing processes (Fig.[5]) with a 
cadence of ^ 10^ yr (Pso,! = 300) for sharp dead zone 
edge models (ARdze < 2AU) and large global viscosity 
(a ^ 0.01). However, the m — 1 mode vortex remains sta- 
bilised for - 5 X lOVr (fto,i = 14000) for thick dead zone 
edge models (2AU ^ ARdze < 5AU) or small global vis- 
cosity (a ^ 0.01). For even thicker dead zone edge models 
(ARdze ^ 5AU), m = 1 mode vortex does not form, only 
m = density perturbation develops. In short, we conclude 
that the horseshoe-shaped vortex (m = 1 mode) formed due 
to the viscosity jump near the disc dead zone edge being at 
^ 50 AU can survive up to the disc life-time (t 5 x 10^ yr), 
if the turbulent viscosity is small (a = 0.001), and the width 
of the dead zone edge is smaller than twice of the local disc 
height (ARdze < 2H). 

To demonstrate the observability of the large vortices, 
we calculated synthetic submillimetre images with a 3D ra- 
diative transfer code using the surface mass density distri- 
bution provided by our hydrodynamic simulations. The dis- 



tribution of the dust, emitting in the submillimetre wave- 
lengths, is assumed to follow that of the gas, owing to the 
dust collecting effect of the pressure bump. The submillime- 
tre intensity distributions given by our ray-tracing models 
- fitting the observed SEDs in the stellar and submillime- 
tre emitting domain (Fig. 11) - show significant azimuthal 
asymmetries (Fig. 10). Comparing our submillimetre predic- 
tions to the observations of three transition discs LkHa333, 
SR21N, and HD 135344B ( [Brown et al.|2009| ), we found as- 
tonishing similarities. Note that there exist further exam- 
ples of discs showing non-axisymmetric submillimetre im- 



ages, such as LkCal 5 (Pietu et al. 2006J, GMAur ( [Hughes 
eraL][2009t, SR24S ([Andrews et al.[[201Q| ), and J1604-2130 



(Williams & Cieza 2011) for which the same phenomenon 



might be responsible. Since the formation process of plan- 
etesimals and planetary embryos could be accelerated by 
the dust collecting and embryo trapping effect of the an- 
ticyclonic vortices caused by the vortical motions and the 
pressure bump fBarge & Sommeria 1995 Klahr & Hen- 



ning 1997), the non-axisymmetric submillimetre images pre- 
sented by Brown et al. ("2009) might be interpreted as snap- 



shots of active planet-forming regions. 

Note that there are other possible routes to generate 
vortices in protoplanetary discs, which also might result in 
m — 1 vortex formation. For the presence of gap-carving 
giant planets the gas is expelled from the planet orbit and 
settles at the gap edges resulting in density maximum which 
is RWI-unstable and develops into vortices as shown by [de| 
jVal-Borro et al.[ ( |2007| ). Another way to generate vortices 



is the baroclinic instability, proposed by ( ^Klahr Boden- 
heimer|2003 ), and recently confirmed by several works (e.g., 



Petersen e t'ar[[2007a|bl [Lesur &: Papaloizou[|20TQl [Lyra k\ 
Klahr 201lf~"^ 

As it was mentioned before, not only enhanced dust co- 
agulation, fragmentation should also happen inside the vor- 
tices if the observed non-axisymmetric features of submil- 
limetre images are indeed due to dust trapped inside them. 
Otherwise, the submilimetre intensity would be weak due to 
the decreasing dust opacity of large grains. Thus, this would 



constitute indirect evidence of elliptic instability ( [Lesur 



Papaloizou 2009) inside vortices leading to 3D turbulence 



in the vortex core, which [Lyra Klahr| ( [2011| ) quantified in 
compressible simulations to be at the level of 10 percent of 
the sound speed. 

Here we have to mention some of the caveats of our 



10 Zs. Regdly^ A. Juhdsz, Zs. Sdndor and C. P. Dullemond 



calculations. Two-dimensional and three-dimensional turbu- 
lence behave quite differently regarding the energy cascade, 
which is form small to large m (direct cascader) in 3D but 
can be form large to small m (inverse cascade) in 2D. 3D 
global disc simulations of Meheut et al. (2010) revealed that 



azimuthal temperature gradients (see, e.g., Petersen et al. 



2007b), but the steep pressure bump and azimuthal den- 



alt hough the RWI results in significant vertical streaming, 
the growth rate and evolution of vortices in 3D is very sim- 
ilar to that of in 2D. This may indicate that the vorticity 
generation by the RWI is strong enough to counter the 3D 
vortex stretching (the mechanism responsible for the direct 
cascade), so that the vortices survive and cascade inversely 
via viscous merging. 

Our simulations cover several million years of the disc 
evolution, yet, the disc photoevaporation acting on the same 
time-scale was neglected for simplicity. Since photoevapora- 
tion process disperse the gas content of the disc, this might 
cause inward migration of the dead zone edge. Thus, photo- 
evaporation process might play an important role in vortex 
evolution, which will be addressed in a forthcoming paper. 



Lin & Papaloizou (2011) have investigated the role of 



self-gravity on vortex instabilities generated at the outer 
edge of the gap opened by a giant planet in locally isother- 
mal disc models. They found that the unstable modes shifts 
to higher azimuthal wavenumbers assuming self-gravity: the 
wavenumber of the most unstable mode increases as Q de- 
creases. For weak self-gravity (Mdisc ^ 0.024 M*), a single 
vortex forms through coalescing, which process is delayed 
as Q decreases. However, for higher disc mass (Mdisc = 
0.031 M*) the final configuration is a vortex pair, m — 2 
mode. Since we used Mdisc = 0.01 M0, our models are in 
the weak self- gravity regime, where the coalescing process 
leads to the formation of m = 1 vortex. 

Self-gravitating disc simulations require the inclusion 
of thermodynamics (the solution of energy equation) too 
{ Pickett et al.|1998| |2000| ^Gammie 2001 ) , otherwise the gas 



concentration would be over predicted. Lyra et al. (2009) 



found that assuming adiabatic equation of state for the gas, 
the pressure bump (that gives rise to RWI) sharpens con- 
siderably compared to their locally isothermal models, due 
to the high temperatures associated with compression. As a 
result, RWI might be triggered easier in adiabatic than in 
isothermal model. Petersen et al. (2007b) investigated the 



baroclinic feedback, in which a vortex enhances azimuthal 
temperature gradients to reinforce the vortex itself, in global 
disc simulations with anelastic-gas model. They concluded 
that this baroclinic feedback mechanism producing stronger 
local temperature gradients results in long-lived vortices if 
the radial temperature gradient is sufficiently large and ei- 
ther the radiative cooling time is small, or the thermal dis- 
sipation is high. On the other hand, vortices in the course of 
vortex coalescing may grow massive enough and reach the 
Jeans length, which provides a natural barrier to growth 
as shown by Mamatsashvili & Rice (2009). In their shear- 



ing sheet simulations, assuming a simple cooling law with a 
constant cooling time and including self-gravity, the vortex 
starts to shock and radiate away the excess vorticity around 
the Jeans length. As a result, the vortices can not grow above 
the Jeans length and short-lived, if there is no continuous 
vorticity production like for RWI, for which case the shear 



is continuously being converted into vorticity (Lyra et al 



sity perturbations formed near the dead zone edge. In short, 
investigating the vortex evolution and longevity in a self- 
gravitating global disc model, considering non-isothermal 
and compressible-gas model, would be important, and will 
be addressed in a forthcoming paper. 

Finally, worth noting that multiple viscosity jumps 
might exist in the disc, e.g., at the inner edge of the dead 
zone or at the snow line dKretke Liii||2007t , which could 
also result in the formation of large-scale vortices. In a non- 
magnetised disc its self-gravity causes outward angular mo- 
mentum transport via quasi-steady self-gravitating turbu- 
lence, which results in the generation of effective viscosity 
(|Pringle|[T98Tt. Since only the outer disc {R >100AU) is 



2009). In our simulations, however, the vorticity is not gen- 



erated by the local advective heat transport enhancing the 



expected to become gravitationally unstable in a standard 
protoplanetary discs, the self- gravity caused effective vis- 
cosity has a gradient there. In such systems the large-scale 
vortex formed beyond ~100 AU might result in vortex aided 
planet formation process at extremely high distances, such 
as observed around Fomalhaut ("Kalas & et al.,"'2008). As a 
conclusion, planets might be assembled in diverse fashion 
by these "planetary factories" located at different distances 
to their host star. In protoplanetary discs where the RWI is 
excited at smaller radii than we investigated, the develop- 
ment of large-scale vortices might lead to the formation of 
a planetary system similar to the Solar System. 
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APPENDIX A: ON THE NUMERICAL 
PRECISION 



As mentioned in Sec. |2.3| the GPU version of FARGO code 
(GFARGO) is only capable of single precision calculation 
at the moment, which obviously results in the degrada- 
tion of precision. In order to be convinced that the results 
of GFARGO (single precision) are appropriate, we ran ad- 
ditional simulations by FARGO (double precision) for the 
model with a 0.1 a^od 0.001, i?dze 50 AU, Ai?dze 
1 AU. In these simulations we applied the same boundary 
conditions, numerical resolutions [Nr — 256, A^(/, = 512), 
and disc masses for GFARGO and FARGO. The surface mass 
density distributions plotted at identical time steps for both 
simulations are shown in Fig. |Al| 

During both simulations, the RWI is excited with m = 6 
azimuthal wavenumber at 3 x 10^ yr (Pso,! = 8). The subse- 
quent vortex coalescing and vortex growth are qualitatively 
equivalent. However, the azimuthal angle of the m — 1 mode 
vortices at 10^ yr (Pso,! = 300) and 1.8 x 10^ yr (Pso,! = 500) 
differs for the single and double precision simulations. This 
is because of the fact that m — 1 mode vortex slowly drifts 
inward (see, e.g., Fig.[4| due to the extra pressure caused by 
the density-enhancement at the dead zone edge, causing pro- 
grade precession of the vortex. Since the FARGO simulation 
results in slightly higher density peak than the GFARGO 
simulation (by ^ 5% higher in average), the vortex drifts 
faster for the double precision, FARGO simulation. Conse- 
quently, the azimuthal angles of the vortices are different. 

In order to investigate the vortex mode oscillation pre- 
sented in Sect. |2.2| we calculated the azimuthal spectral 
power of the surface mass density distribution against time 



(Fig. A2 ). As one can see, the evolution of the spectral power 



in the ^ m ^ 6 modes are very similar for the two sim- 
ulations. Note, however, that although the first and second 
vortex-decay occur at the same time for both simulations 
(Pso,! — 600 and Pso,! — 1000, respectively), the time re- 
quired to reform the stable m = 1 mode is slightly longer 
for the double precision simulation. 

In short, we conclude that the lack of double precision 
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Figure Al. Surface mass density evolution in a fast evolving 
model (i.e. model with a = 0.1, oimod = 0.001, Rdze = 50 AU, 
Ai^dze = lAU). The simulations were done by single precision 
GFARGO code {left column) and double precision FARGO code 
(right column) with exactly the same boundary conditions and 
numerical resolution. The surface mass density (shown with verti- 
cal colour bars for each frame) is measured in g cm~^ . The frames 
were saved at exactly the same time steps for both simulations. 
P50,i means the number of orbits at 50 AU assuming a Solar mass 
central star. 



capability of GFARGO is not dangerous on a million year 
time-scale in our simulations. 



APPENDIX B: ON THE EXCITATION OF RWI 

In order to make sure that no important feature is being 
missed during the excitation of the RWI, we ran several mod- 
els where the RWI was seeded by conventional methods. To 
provide seeds for the RWI, we added low level perturbations 




400 600 800 1000 

Time (Pso.i) 




Figure A2. Azimuthal spectral power of the surface mass density 
distribution against time measured in Pso,! for simulations shown 
in Fig.[Al] Results given by GFARGO and FARGO are shown in 
the upper and lower panels, respectively. The Fourier amplitudes 
are normalised by the total power measured in the m = — 6 
modes. 



to the surface mass density when the threshold for the RWI 
had been reached (i.e., when significant density-bump had 
been developed at the dead zone edge) . In the first series of 
simulations, the perturbation was calculated as a superposi- 



tion of sines, similarly to Eq. (10) presented in ,Meheut et al. 
(|20T0|): 



S(i?,( 



S(i?,(/)) + |S(i?,(/))|esin(27ri?) x 

[sin{(j)) + sin(2(/)) + sin(3(/)) + sm(5(/))] , (Bl) 



where S(i^, 0), and S(i?, 0) are the unperturbed, and per- 
turbed surface mass density, respectively. The magnitude of 
the perturbation is set by e. In the second series, the surface 
mass density were perturbed by normal distribution white 
noise. While the former type of perturbation has a drawback 
that it excites only the low-m modes as does a small mass 
orbiting planet, large part of the spectrum is being excited 
by the latter type of perturbation. We also tested the RWI 
excitation by small mass planet orbiting at 50 AU. 

The RWI excitation was investigated in model where 
parameters a — 0.01, amod — 0.01, Rdze = 50 AU, and 
AR = 1 AU were chosen, while the grid resolution was set to 
Nr = 256 and =512. The magnitude of perturbations 
were in the order of e = 10~^. The simulations were done by 
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FARGO, because it has smaller numerical noise in the radial 
velocity component than in GFARGO simulations. 

First, the threshold for the RWI was let to develop, 
i.e. a significant radial density gradient with magnitude of 
AE/E ^ 2 at the dead zone edge. Then one of the above 
described perturbation was applied. To ascertain whether 
the RWI is not excited by the numerical noise of FARGO, 
we ran models without perturbation. In this unperturbed 
model, we found no RWI excitation. 

Figures |B1| shows the evolution of azimuthal spectral 
power of the surface mass density followed as long as m = 1 
mode vortex had been fully fledged (by Pso,! = 150) for 
models using different type of perturbations: superposition 
of sines in surface mass density with magnitude of e = lO"'^ 
{upper panel)] normal distribution white noise in surface 
mass density with magnitude of e = 10~^ {middle panel)] 
small mass (Mpi = 10~^M*) planet orbiting at 50 AU {lower 
panel). As one can see the evolutions of spectral power are 
very similar for any type of perturbations. Independent of 
the perturbation type applied, we found that i.) the RWI is 
always triggered after Pso,! — 5 — 10 orbits at 50 AU; ii.) the 
most unstable mode of RWI is m = 5 or m = 6 at the onset 
of vortex formation; iii.) following the evolution of the vor- 
tices, m — 1 mode (horseshoe-shaped) vortex always forms. 
Consequently, the mode of the excitation of RWI has no im- 
portance regarding our main findings that long-lived m = 1 
mode vortex appears near the disc dead zone edge. 




Figure Bl. Evolution of the azimuthal spectral power of surface 
mass density distribution for models ol — 0.01, oij^iod 

= 0.01, 

Ai^dze = 1 AU followed as long as m = 1 mode vortex had 
been fully fledged (at P5o,i = 200). The RWI was seeded by su- 
perposition of sines {upper panel), by normal distribution white 
noise {middle panel) in surface mass density, and by small mass 
(Mpi = 10-'^ M*) planet orbiting at 50 AU (/ ower vanel). The 
Fourier amplitudes are normalised by the total power measured 
in the m = — 6 modes. 



